library(aqp)
library(soilDB)
library(compositions)
library(reshape2)
library(latticeExtra)
library(plyr)

# get sand, silt, clay data from the named series
# return weighted mean of top 10cm
get.data <- function(series) {
  x <- fetchKSSL(series)
  s <- slab(x, pedon_key ~ clay + sand + silt, slab.structure=c(0,10), strict = FALSE, slab.fun = mean, na.rm=TRUE)
  s <- dcast(s, pedon_key ~ variable, value.var = 'value')
  ssc <- na.omit(horizons(x)[, c('sand', 'silt', 'clay')])
  return(ssc)
}

# combine original texture + simulated textures from normal composition
sim.data <- function(ssc, n.sim=1000) {
  # convert to compositional class, note range is now [0,1]
  ssc.acomp <- acomp(ssc)
  # simulate normally-distributed composition based on data
  ssc.sim <- rnorm.acomp(n=n.sim, mean=meanCol(ssc.acomp), var=var(ssc.acomp))
  # convert to data.frame and original range of [0,100]
  ssc.sim <- as.data.frame(unclass(ssc.sim) * 100)
  # return stacked original + simulated data
  return(make.groups(original=ssc, simulated=ssc.sim))
}


# get data for named series
s.list <- c('clarksville', 'vista', 'auburn', 'cecil', 'drummer', 'capay')
d <- lapply(s.list, get.data)
names(d) <- s.list
d <- ldply(d)

# simulate from a normal composition, using mean and var-cov matrix from original data
d.sim <- ddply(d, '.id', function(i) sim.data(i[, -1]))
m <- melt(d.sim, id.vars=c('.id', 'which'))
# plotting style
cols <- brewer.pal(n=3, 'Set1')
tps <- list(superpose.line=list(col=cols), superpose.symbol=list(col=cols, cex=0.5, pch='|'))

# plot: panel according to original | simulated 
p.1 <- densityplot(~ value | which + .id, groups=variable, data=m, xlab='Sand, Silt, Clay (%)', auto.key=list(columns=3), par.settings=tps, scales=list(alternating=3), panel=function(...) {
  panel.grid(-1, -1)
  panel.densityplot(...)
})
p.1 <- useOuterStrips(p.1)

# plot: panel according to texture component
p.2 <- densityplot(~ value | variable + .id, groups=which, data=m, xlab='Sand, Silt, Clay (%)', auto.key=list(columns=2), par.settings=tps, scales=list(alternating=3), panel=function(...) {
  panel.grid(-1, -1)
  panel.densityplot(...)
})
p.2 <- useOuterStrips(p.2)

# display plots
print(p.1)

print(p.2)

## viz on texture triangle
texture.triangle.low.rv.high(d[d$.id == 'clarksville', -1], sim=TRUE, cex=0.5); title('Clarksville')

texture.triangle.low.rv.high(d[d$.id == 'drummer', -1], sim=TRUE, cex=0.5); title('Drummer')

texture.triangle.low.rv.high(d[d$.id == 'auburn', -1], sim=TRUE, cex=0.5); title('Auburn')

texture.triangle.low.rv.high(d[d$.id == 'cecil', -1], sim=TRUE, cex=0.5); title('Cecil')

texture.triangle.low.rv.high(d[d$.id == 'vista', -1], sim=TRUE, cex=0.5); title('Vista')

texture.triangle.low.rv.high(d[d$.id == 'capay', -1], sim=TRUE, cex=0.5); title('Capay')

This document is based on aqp version 1.10-3.